The formula on which the implementation of Reynolds stress in coding is based on is
\[ \begin{aligned} \langle U_i U_j \rangle &= \langle ( \langle U_i \rangle + u_i ) ( \langle U_j \rangle + u_j ) \rangle \\ &= \langle \langle U_i \rangle \langle U_j \rangle + u_i \langle U_j \rangle + \langle U_i \rangle u_j + u_i u_j \rangle \\ &= \langle U_i \rangle \langle U_j \rangle + \langle u_i u_j \rangle \end{aligned} \]
Therefore, the implementation of Reynolds stress is similar to the time averaging of the velocity. \(\langle U_i U_j \rangle\) is the new target. The implementation would be
\[ \langle u_i u_j \rangle_{n} = \frac{(n-1) \left (\langle U_i \rangle_{n-1} \langle U_j \rangle_{n-1}+\langle u_i u_j \rangle_{n-1} \right )+\left ( U_i U_j \right )_{n}}{n} - \langle U_i \rangle_{n} \langle U_j \rangle_{n} \]